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ABSTRACT 

Aims. A numerically stable and accurate advection scheme for modeling shock fronts moving at ultra-relativistic speeds is funda- 
mentally important for understanding the thermodynamics of jets emanating from the vicinity of relativistic objects and the origin of 
the high gamma-rays. 

Methods. In this paper we present a spatially third-order accurate advection scheme that is capable of capturing shock-fronts with 
diverse Lorentz factors. The consistency and accuracy of the scheme are investigated using the internal and total energy formulation 
in general relativity. 

Results. Using the total energy formulation, the scheme is found to be viable for modeling moving shocks at moderate Lorentz fac- 
tors, though with relatively small Courant numbers. In the limit of high Lorentz factors, the internal energy formulation in combination 
with a fine-tuned artificial viscosity is much more robust and efficient. 

We confirm our conclusions by performing test calculations and compare the results with the analytical solutions of the relativistic 
shock tube problem. 
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1. Introduction 

Plasmas moving at relativistic speeds have been observed in 
diverse astrophysical phenomena, such as in supernova explo- 
sions, in jets emanating from around neutron s tars, microquasars 
and in active galactic n uclei (see iLivio L 120041 : iMarscherl 
l2006HFender" et al. L I2OO7I and the references therein). The bulk 
Lorentz factor, F, of the observed jet-plasmas, in most cases, has 
been found to increase with the mass of the central relativistic 
object. Ejected plasmas from around neutron stars have F^s - 
0(1- 3), in microquasars F^qso -0(1-4) and in AGNs/QSOs 
^ACN - 0(5 - 10). In gamma-ray bursts however, plasmas are 
considered to move with a L orentz factor of the order of several 
hundreds (e.g.. |piranLl2006l) . i.e., Tcrb = 0(lOO - 1000). 
Jet-plasmas considered to decelerate via self-interaction or in- 
teraction with the surrounding media in the form of shocks and 
eventually become efficient sources for the production of the ob- 
served high energy gamm a-rays and probably for the energetic 
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cosmic rays ( Sikorairi994l) . 
Modeling the formation of relativistic shock by means of 
HD and MHD simulation has been the of fo cus of numer- 
ous studies during the t he last two de cades (Hawley et tdj, 
■l984b; Alav et a l.l [1991 iFontetal 
20031: IDe Vilhers. HawlevL 12003 



12000; Marti and Mulle 
iGammie et al.L I20()3l 
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20071: iMignone et al.L I2OO7I: iHuieirat et all I2OO8I see also the 
references therein). 

These studies have enriched the field of computational astro- 
physics with numerical techniques and useful strategies for ac- 
curate capturing relativistic shock fr onts using both the internal 
and total energy for mulation (e.g., iDe Vilhers. Hawlevl 120031 : 
iMignone et al.Ll200 7). Nevertheless, most of the methods appear 
to encounter severe numerical difficulties whenever the Lorentz 
factor becomes large, i.e., F > 5. The reason is that small veloc- 
ity en^ors are magnified by F^ as they enter the other equations, 
hence inducing super-proportional errors in the estimation of the 
other variables. Recipes such as further reducing the grid spac- 
ing and accordingly the time step size my lead to stagnation of 
the solution procedure, particularly if the methods used are of 
time-explicit type. 

The robustness of time-implicit methods at capturing shocks 
propagating with high Lorentz factors has been barely inves- 
tigated, hence the purpose of the present paper. Moreover, we 
show that our time-implicit method in combination with a mod- 
ified third order MUSCL advection scheme is capable of model- 
ing shocks propagating with Lorentz factors F > 10 with time 
step sizes corresponding to Courant numbers larger than 1 . 
In Sec. 2 we describe the hydrodynamical equations solved 
in the present study. The solution method relies on using the 
3D axi-symmetric general relativistic implicit radiative MHD 
solver (henceforth GR-I-RMHD), which has been de s cribed 
in details in a series of ar ticles (e.g., iHuieirat et al.L 120081: 
iHuieirat & Thielemann.ll200 9). In Sec. (3) we just focus on sev- 
eral new aspects of the advection scheme. The results of the ver- 
ification tests is presented in Sec. (4) and end up with Sec. (5), 
where the results of the present study is summarized. 
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2. The 1D general relativistic hydrodynamical 
equations 

In the present study we consider the set of Euler equations in 
one-dimension and in flat spacetime: 



dtD + V,. ■ (D VO = 
dtMr + V, ■ (MrV) = -drP 
d,&°' + Vr-([S'°' + P]V) ^0, 



(1) 



where Vj - = 



g, D(= pu') is the relativistic density, 



Ml is the radial momentum: Mr = Dut, where D = Dh, u' is the 
time-like velocitjQ , V - u'^'/u' is the transport velocity, "h" is 
the relativistic enthalpy h = c^ -i- e H- P/p. = (u')^ph - P is 
the total energy which is the sum of kinetic and thermal energies 
of the gas. "P" is the pressure of ideal gas: P = (y - l)ps, s and 
7 denote the internal energy and adiabatic index, respectively. 

The reader is referred to Sec. (2) of Hujeirat et al. (2008), 
where the general relativistic hydrodynamical equations and 
their derivations in the Boyer-Lindquist coordinates are de- 
scribed. However, we continue to use these coordinates, though 
in the limit of flat space. 

In the case that the mechanical energy is conserved, then the 
total energy reduces to an evolutionary equation for the internal 
energy: 



cd 

5t£'* + V,. ■ (£'^¥0 = -(r - 1) — [5,u' + Vr ■ (u' Y)] 



(2) 



where = u'P/(y - 1). 

Furthermore, there is a remarkable similarity of how Eq. ^ 
deals with S'^ and u'. Taking this similarity into account, the 
equation can be re-formulated to have the following form: 



8,6'^ + Vr ■ (£''V'") = -rD(V,- ■ V"), 



(3) 



where = Dlog[£'^(u'F ']. 

In Equation ([T]i, the energy equation describes the time-evolution 
of the total energy fi'°*. Our numerical solution method relies 
first on selecting the primary variables and calculating their cor- 
responding entries in the Jacobian. We then iterate to recover 
the correct contributions of the depending variables, such as the 
transport velocity and pressure. Therefore, the advection opera- 
tor can be decomposed into two parts as follows: 



5t£'°' -t- V,- ■ (fi*°'V') = Vr ■ (P V). 



(4) 



This equation is solved as follows: Solve for £'°' using the pres- 
sure from the last iteration level. Once and "D" are known, 
we may proceed to update the pressure as follows: 



P = 



£'°' - u'D 

y-1 V / 



(5) 



This value of "P" is then inserted again into the RHS of Equation 
(|4|i to compute a corrected value for 

The effect of iteration can be reduced on the cost of using pres- 
sure values from the last time step using the following alternative 
formulation: 

5tfi'°' -I- V,. ■ (fi'°'V') = 5tP, (6) 

' u' is also the general relativistic Lorentz factor, which reduces to T 
in flat space. 



where ^ (u'fph. 

Thus, once is found, the pressure can be computed from the 
known conservative variables as follows: 



P = 



u*D 



r 



(uO 



t^2 



(7) 



On the other hand, using the internal energy formulation for 
modeling moving shocks requires the inclusion of an artificial 
viscosity for calculating their fronts accurately. This modifica- 
tion is necessary in order to recover the loss of heat generally 
produced through the conversion of kinetic energy into internal 
energy at the shock fronts. 

In this case, the RHS of Equation (|2]i should be modified to in- 
clude the artificial heating term: 



(8) 



where rjart is the artificial viscosity coefficient. 
The inclusion of 2a,/ implies an enhancement of the effective 
pressure. Therefore, the thermodynamical pressure, and also the 
enthalpy, should be modified to include an artificial pressure of 
the form: 



Ptot = P + Pan = P + 77art5;.(u'VO. 



(9) 



Note that the artificial pressure enters the momentum equation 
in the form of VP„r(, which scales as ~ A(riar,V). This is equiva- 
lent to activating a second order viscous operator at the shock 
fronts, whose effect is then to transport information from the 
downstream to the upstream regions, so to enhance the stabil- 
ity of the transport scheme in such critical transitions. 



2. 1 . Determining the Lorentz factor 

Our test calculations showed that the Lorentz factor can be best 
determined from the normalization condition, in which both the 
conservative variables and transport velocities are involved. To 
clarify this point, the normaUzation condition reads: 



-1 



A/fl ..lip ,, 



§[M, + VMr + V^Me + V^M^] 

D g" 



(10) 



where D - phu' and which can be directly determined from the 
continuity and the energy equations, is the transport velocity 
and g''^ are the elements of the metric. 

The final form of Eq. ( fTOt is the quadratic equation: a (u')^ -i- 
b u' + c - 0, from which the Lorentz factor can be determined 
completely. 

We note that since the state of the conservative variables 
D, M^,fi'' depends strongly on the transport velocity V, it is 
therefore necessary to consider V when computing the Lorentz 
factor u' as described in Eq. ( fTOl l. 



3. A modified third order advection scheme for high 
Lorentz factors 

Most time-implicit advection schemes generally rely on upwind- 
ing to stabilize the transport near critical fronts. Transport of 
quantities here can be mediated with CFL numbers larger than 
unity. Therefore, the advection scheme employed should be in- 
dependent of the time step size. Indeed, the monotone upstream 
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centered schemes for conservation l aws, or simpl y the MUSCL 
schemes, obey these conditions (see' Hirschlll98 8', and the refer- 
ences therein). However, in order to model moving shocks with 
high Lorentz factors accurately, it is necessary to modify the 
scheme, so that the information-flow from the downstream to 
upstream regions should be further limited in accordance with 
the normalization condition, i.e., causality condition. 
A classical advection scheme of the MUSCL-type gives the fol- 
lowing interface values: 



,SR,SL 



SR 
4j-l/2 

SL 
4j-l/2 



= qj_i - + ^)Aqj_i + (1 - ^)Aqj_2] : if Vj < 
= qj + - ^)Aqj + (1 + /f)Aqj+i] : if Vj > 0, 



(11) 



where "q" is the transported physical quantity, /c is a switch 
offyon parameter that specify the accuracy needed and Aqj - 
qj - qj+i . Note that the scheme is of second order for k - - 1 , 0, 1 
and of third order for k - 1/3. 

In order to enable an accurate capturing of shock fronts propa- 
gating with high Lorentz factors on a non-uniform grid distribu- 
tion, the following modifications have been performed: 



if Vj<Q: 



•cr)Aqj + (1 

^RR\2 



if 



,j qj-i - fL(l -cr)Aqj + (1 +cr)Aqj_i], 
where o-^2k {Ac^^ ■ A<^^)l[{Ac^^f + (Aqj'^)^ -i- e]. 
Vj>Q: qfL = qj + i[(l + (r)Aqj + (1 - (r)Aqj+i], 



(12) 



where (t^Ik (Aq]'^ ■ Aq[^)/[(Aq]'^)^ -i- (Aq[^)^ + e]. 



where 



Aqj^R = (xP, - x™)(qj_2 - qj-i)/(x]"2 • 

Aqf = (x-2-Ti)(qj-i 

xr.i)(qj-i-qj)/(xf_i 

(xPj - xp(qj - qj+i)/(x]" - 



Aq°L = (x| 
Aqj-L 



-xp 



(13) 



and where e is a small number set to avoid division by zero. 



where ^ is an additional weighting function and Lq. 
LUIS corrections of "q" at the interface rj (see Fig. |2] which 
are computed in the following manner: 



Lq 



SR 



l=j-2 41 lli^k,k=j-2 r"'-i 



k=j 



if Vj<0 



SI v-'l=i+2 -i-rk=i+2 



if 



Vj > 0. 



(15) 



However, we note that high order advection schemes gener- 
ally degenerate into lower ones at the shock fronts or discon- 
tinuities, in order to maintain monotonicity of the advection 
scheme, which in turn enhances the production of numerical en- 
tropy/diffusion at these critical transitions. Therefore, in most 
of the test calculations presented here it was necessary to have 
^ > 1 /2, or even very close to unity to enable convergence. 

4. Results 

In this section we present results of various verification tests us- 
ing different energy formulations in combination with the advec- 
tion scheme described in the previous section. The widely used 
hydrodynamical test is the one-dimensional relativistic shock 
tube problem (henceforth RSTP). The advantages of this test 
problem is that the obtained numerical solutions can be com- 
pared with the corresponding analytical solution with arbitrary 
large Lorentz factors. 

For further details on the background of this test problem and 
the way the analytical solutions are obtained we refer the reader 
to divlarti and Muller, 2003). 

The RSTP is an initial value problem, in which the solution de- 
pends strongly on the initial conditions. So the domain of cal- 
culations is divided into a right and left regions (see Fig.|2]l. In 
the left region < x < Xc : we set Pl = 40/3, pl - 10 and in 
the right region Xc < x : we set Pr = (2/3) x 10^^, pr = 1. The 
initial velocity is set to vanish everywhere. 



1, ; ?J 



j-1 



r. 



j+2 



^j+1 



J 



J- 



Fig. 1. A schematic description of several finite volume cells, their 
boundaries and the location of the scalar quantities "q" and the veloci- 
ties. 



dx 



Rr 



Fig. 2. Grid spacing versus position. The actual changes of the vari- 
ables occur in the interval [x^, < x < Xi,, where the grid spacing is 
uniform. The initial conditions read: for x < x^, : Pl = 40/3, pl = 10 
and for x > x,. : Pr = (2/3) x 10"* and p« = 1. 



The accuracy of above modified scheme may be further in- 
creased by incorporating the Lagrangian Upwind Interpolation 
Scheme (LUIS). The LUIS strategy is based on by constructing 
a Lagrangian polynomial of third or fourth order whose main 
weight is shifted to the right or left, depending on the upwind 
direction. The final combined LUIS-MUSCL scheme reads as 
follows: 

if y, < : q^^ = ^q^'* -h (1 - ^Lq^"^, 

if Vj>0: qf^ = ^qp + (I - ^)Lqp, ^ ' 



Using the total energy formulation (TEF) we show in Figures 
(l3lfT3]) the profiles of the velocity and the coiTesponding Lorentz 
factor for different distributions of the initial conditions. 
In Fig|3]we show the profiles of the velocity and Lorentz fac- 
tor after 0.2 sec using the total energy formulation (see Eq. [U. 
The interval [0.35 < x < 0.75] has be divided into 400 equally 
spaced cells, whereas 220 cells where used to cover the exter- 
nal intervals, where the variables continue to acquire their initial 
values. The advection scheme described in (Eq. [T4l i has been 
employed to achieve 3'^'' order spatial accuracy and the damped 
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Crank-Nicols on method is used to achiev e second order tempo- 
ral accuracy ("Hujeirat. Rannacheri 1200 Ih . The time step size is 
set to have a maximum that corresponds to CFL= 0.4. 
To test the robustness of the solution procedure, we have in- 
creased the Pl by a factor of 100, while the other variables 
remained unchanged. This increase of the pressure yields the 
Lorentz factor: u' » 3.6 (Fig. |4|i. However, to enforce conver- 
gence, it was necessary to increase the number of iteration per 
time step considerably. 

Comparing the results with those of Fig. (|5]l, we see that the 
solution obtained is sufficiently accurate, so that doubling the 
number of grid points and decreasing the time step size by half 
did not improve the accuracy of the scheme significantly. 
From Figs. (|6]l and (|7J, we see that incorporating artificial vis- 
cosity and artificial heating while using the TEF has a similar 
effect as that of the normal numerical diffusion. In both cases 
the velocity at the shock front attains lower values than expected 
from the analytical solution. 

We note that within the context of time-implicit solution meth- 
ods, the total energy method appear to have a limited capability 
to capture shocks moving with high Lorentz factors (Fig.lSJ. The 
reason is that, unlike the internal energy Ein, the kinetic energy 
Ekin increases with Lorentz factor as F-. As a consequence, the 
convergence rate decreases with Ein/Ekin- In this case it is nec- 
essary to decrease both the time step size and the grid spacing 
in addition to increasing the number of iteration per time step. 
This procedure may easily lead to a stagnation of the solution 
method. 

On the other hand, our time-implicit solution method appears 
to be more robust if the internal energy formulation is used. In 
Figures (l9lfT3ll. the results of several test calculations are shown, 
which agree well with the analytical solutions. Moreover, the 
IFF is sufficiently stable, so that running the calculations with 
CFL > 1 does not appear to hinder the convergence of the im- 
plicit solution procedure (see Figs. [T2] and [T3Tl. 
The major draw back of the IFF is the necessity to fine-tune 
the coefficients of the artificial viscosity in order to obtain the 
correct Lorentz factor that fits well with the analytical solution. 
Moreover, these fine-tuning procedures of the parameters must 
be repeated each time the ratio of the pressure Pl/Pr is changed. 
This is a consequence of the conversion of kinetic energy into 
heat at the shock fronts which cannot be determined a priori and 
uniquely by solving the internal energy alone. 
Nevertheless, using the IFF, non-linear physical processes can be 
easily incorporated into the implicit solution procedure. Unlike 
in the TEF case, such processes are direct functions of the inter- 
nal energy and therefore it is much easier to compute their entries 
and incorporate them into the corresponding Jacobian. This has 
the consequence that calculations can be performed also with 
CFL > 1, while still maintaining consistency with the original 
physical problem. 

Finally, in Fig. flAi we show the profiles of V and u' using the 
compact internal energy formulation (CIEF, Eq.Q with high ac- 
curacies. The method is found to display overshooting and un- 
dershooting both at the contact discontinuity and the shock front, 
whose appearing was difficult to prevent through enlarging the 
number of grid points, decreasing the time step, enhancing the 
number of iteration per time step or even by activating the artifi- 
cial viscosity. 

It should be noted here that since in the limit of large Lorentz 
factors, the CFL number depends weakly on the sound speed 
(see Fig. [TSl l, increasing the global iteration per time step has 
more significant effect on convergence than merely decreasing 
the time step size. 




Fig. 3. The profiles of the velocity V (left) and the Lorentz factor 
u' (right, dashed line connecting squares) overplotted on the corre- 
sponding analytical solutions (continuous blue lines) after t = 0.2. 
The numerical solutions have been obtained using the total energy for- 
mulation (TEF), the initial pressure (Plq = 40/3, Pl = 10), (Pr = 
(2/3) X 10"*, PR = 1), 700 non-uniformly distributed grid points and 
CFL=0.4. The advection scheme employed here is of third order spatial 
and second order temporal accuracy. 




Fig. 4. As in Fig (3), the profiles of V and u' have been obtained us- 
ing the TEF, an initial pressure Pl = 10~Plo and CFL= 0.4. 620 non- 
uniformly distributed grid points have been used 




Fig. 5. As in Fig (3), the profiles of V and u' have been obtained using 
the TEF, an initial pressure Pl = IO^Plo, CFL= 0.4 and 1300 non- 
uniformly distributed grid points. 



5. Summary and conclusions 

In this paper we have presented the results of different numer- 
ical studies for modeling the propagation of ultra-relativistic 
shocks using the internal energy (lEF), the total energy (TEF) 
and pseudo-total energy (PTEF) formulations. 
Our test calculations show that TEF is most suited for modeling 
the propagation of shocks with low to moderate Lorentz factors, 
though with low CFL numbers, hence more suitable for time- 
explicit solution methods. 

On the other hand, although the PTEF is similar to the TEF, it 
was found that a larger number of iterations per time step was 
required to maintain convergence. 

Using the lEF however, it was verified that the implicit solution 
procedure converges relatively fast even for a Lorentz factor 
u' » 1 and a CFL number of CFL > 1 . 
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Fig. 6. As in Fig (3), the profiles of V and u' have been obtained us- 
ing the TEF, an initial pressure Pl = IO^Plo, a CFL= 0.4. An upwind 
advection scheme of first order spatial accuracy has been used. 



Fig. 10. Similar to Fig (3), the profiles of V and u' have been obtained 
using the lEF with the initial left pressure Pl = IOPlo- The rest of 
variables remained as in Fig. (1) unchanged. 



— Analytic solui 
.- -. 620 grid poii 







Fig. 7. As in Fig (3), the profiles of V and u' have been obtained us- 
ing the TEF with a heating source due to artificial viscosity, an initial 
pressure Pl = 10"Plo, CFL= 0.4 and a spatially third order accurate 
advection scheme. 



— Analytic solui 
■- -■ 700 grid poii 





Fig. 8. As in Fig (3), the profiles of V and u' have been obtained us- 
ing the pseudo-total energy formulation (PTEF, see Eq. [6](, an initial 
pressure Pl = IO^'Plo and CFL= 0.4. 




ra 1.25 ■ 

£ 1.20 . 




Fig. 9. As in Fig (3): the profiles of V and u' have been obtained using 
the internal energy formulation (IFF, see Eq.O, an initial pressure Pl = 
Plo and CFL= 0.4. 



Fig. 11. As in Fig (3), the profiles of V and u' have been obtained using 
the IFF, an initial pressure Pl = IO^Plo and aCFL =1.5 after r = 0.215. 





Fig. 12. As in Fig (3), the profiles of V and u' have been obtained using 
the IFF, an initial pressure Pl = IO'Plo and a CFL= 2. 



— Analytic solutii 
.- -. 700 grid point; 





Fig. 13. As in Fig (3), the profiles of V and u' have been obtained using 
the lEF, an initial pressure Pl = IO^'Plo and a CFL= 4. 



into the Jacobian. 



It should be stressed here that the lEF has a major draw back, 
as it requires a fine-tuned parameters of the artificial viscosity 
to obtain the correct Lorentz factor at the shock fronts. On the 
other hand, when modeling astrophysical plasmas, the lEF is 
favorable over the TEF, as cooling and heating processes are in 
general direct functions of the internal energy, thus it is easier to 
construct the corresponding coefficients and incorporate them 



Furthermore, turbulence is an important property of most 
astrophysical fluid flows. Thus the artificial viscosity in such 
flows naturally will be overwhelmed by turbulent diffusion, 
hence relaxing the fine-tuning problem of the artificial viscosity. 

Irrespective of the method used for solving the energy equation, 
the Lorentz factor is found to be best computed from the 
normalization condition in which the transport velocities and 
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Fig. 14. As in Fig (3), the profiles of V and u' after t = 0. 175 have been 
obtained using the compact internal energy formulation (CIEF, see Eq. 
O and the initial pressure Pl = IOPlo- 



10"^ U — . — — . — — . — — . — — . — 

10 10 10 10 10 10 

Fig. 15. fH = Vs / c for the Newtonian case (Vs = yP/p; dashed line) and 
for the relativistic case (Vs = yP/ph; solid line) versus temperature. 

momenta are used. 

We have also presented a modified MUSCL advection scheme 
of third order spatial accuracy, which has been constructed to 
enable accurate capturing of relativistically moving shocks. The 
accuracy of the scheme can be further enhanced by incorporat- 
ing the Lagrangian Upwind Interpolation Scheme (LUIS). 
On the other hand, using this advection scheme for modeling 
astrophysical fluid flows may unnecessary increase the compu- 
tational costs, in particular when 3D axi-symmetric calculations 
with high resolution are required. 

Noting that the appropriate parameters of the artificial viscosity 
cannot be determined a priori, obtaining the correct Lorentz fac- 
tor may turn out to be a difficult task. In the case of stan ding 
shocks, the hierarchical solution scenario jHuieiratl l2005h can 
be employed to gradual enhance the accuracies of the scheme to 
obtain the correct Lorentz factors at the shock fronts. 
In a forthcoming paper, we intend to study the formation of 
strong shocks in curved spacetime near the surface of ultra- 
compact neutron stars. 
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